﻿// ******************************************************************************************************
// <copyright file="Kriging.cs" company="Caroso.inc">
//     Copyright (c) Carlos Osorio All rights reserved.
// </copyright>
// ******************************************************************************************************
// The contents of this file are subject to the Mozilla Public License Version 1.1 (the "License"); 
// you may not use this file except in compliance with the License. You may obtain a copy of the License at 
// http:// www.mozilla.org/MPL/ 
// Software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTY OF 
// ANY KIND, either express or implied. See the License for the specificlanguage governing rights and 
// limitations under the License. 
// 
// The Original Code is this .mapWindow 
// 
// The Initial Developer of this version of the Original Code is Carlos Osorio carosocali@gmail.com
// 
// Contributor(s): (Open source contributors should list themselves and their modifications here). 
// Change Log: 
// Date           Changed By      Notes
// 22 Nov 2010    Carlos Osorio   Inital upload t
// ******************************************************************************************************
namespace GeostatisticalTool.Classes
{
    using System;
    using System.Collections;
    using DotSpatial.Data;
    using DotSpatial.Topology;
    using DotSpatial.Projections;
    using DotSpatial.Controls;
    using System.Collections.Generic;
    using System.Threading.Tasks;

    /// <summary>
    /// Main class  for Kriging Interporlation
    /// </summary> 
    public class Kriging : AbstractInterpolator, IKriging
    {
        #region propierties
        /// <summary> Number of points</summary> 
        //  private int npoints = 0;

        /// <summary>This is the Binding</summary> 
        private double[,] mlags;

        /// <summary>Matrix with the coordinates of all points</summary> 
        // private double[,] cpoints;

        /// <summary> Distance max </summary> 
        private double maux = 0;

        /// <summary>Max value of Semivariance</summary> 
        private double mauxV = 0;

        /// <summary>Initial lag.</summary> 
        private double lag = 0;

        /// <summary>Model defintion</summary> 
        private IKrigingModel mmodel;

        private double[] vectorG;
        private double[] weights;



        #endregion

        /// <summary>
        /// Initializes a new instance of the Kriging class
        /// </summary>
        /// <param name="shapeLayer">Shapefile type point</param>
        /// <param name="idField">Position of field that contains the Z value</param>
        public Kriging(IFeatureSet shapeLayer, int idField) :
            base(shapeLayer, idField)
        {
            this.method = "Kriging";
            var modelLineal = new ModelLineal
                                  {
                                      C0 = 0,
                                      C1 = 1,
                                      Range = 1
                                  };
            // this.CreatePoints(shapeLayer, idField);
            this.CalculateDistances();
        }




        /// <summary>
        /// Initializes a new instance of the Kriging class
        /// </summary>
        /// <param name="shapeLayer">Shapefile (POINTS)</param>
        /// <param name="idField">Position of field that contains the Z value</param>
        /// <param name="model">This class provides the model that will be used</param>
        public Kriging(IFeatureSet shapeLayer, int idField, IKrigingModel model)
            : base(shapeLayer, idField)
        {
            this.method = "Kriging";
            if (model != null)
            {
                this.mmodel = model;
            }
            this.CalculateDistances();
            //  this.CreatePoints(shapeLayer, idField);
            ////  GenerateMatrix(ref msemiVarTeo, ref inverMatrix, ref mdistances);
        }

        #region fields


        public IKrigingModel KrigingModel
        {
            get { return mmodel; }
            set { mmodel = value; }
        }


        /// <summary>Gets a value indicating max Distance</summary> 
        public double MaxDistance
        {
            get { return this.maux; }
        }

        /// <summary>Gets a value indicating max variance</summary> 
        public double MaxVariance
        {
            get { return this.mauxV; }
        }



        /// <summary>Gets or sets the  length of the lag</summary> 
        public double ValueLag
        {
            get { return this.lag; }
            set { this.lag = value; }
        }

        #endregion



        /// <summary>
        /// Calculate the distastance between two points
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>Return the distance</returns>
        public double Mdistances(
            int i, int j)
        {
            return Math.Sqrt(((this.Kpoints[i].X - this.Kpoints[j].X) * (this.Kpoints[i].X - this.Kpoints[j].X)) +
                ((this.Kpoints[i].Y - this.Kpoints[j].Y) * (this.Kpoints[i].Y - this.Kpoints[j].Y)) + ((this.Kpoints[i].Z - this.Kpoints[j].Z) * (this.Kpoints[i].Z - this.Kpoints[j].Z)));

            //    return Math.Sqrt(((this.cpoints[i, 0] - this.cpoints[j, 0]) * (this.cpoints[i, 0] - this.cpoints[j, 0])) + ((this.cpoints[i, 1] - this.cpoints[j, 1]) * (this.cpoints[i, 1] - this.cpoints[j, 1])));
        }

        /// <summary>
        /// Calculate Semivariances
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>the semivariance</returns>
        public double MsemiVarEmp(
             int i, int j)
        {
            return 0.5 * ((this.Kpoints[i].W - this.Kpoints[j].W) * (this.Kpoints[i].W - this.Kpoints[j].W));
            //   return 0.5 * ((this.cpoints[i, 2] - this.cpoints[j, 2]) * (this.cpoints[i, 2] - this.cpoints[j, 2]));
        }



        /// <summary>
        /// CalculateLags
        /// This method calculates the semivariance taking into account the lag
        /// </summary>
        /// <param name="lag">Distance. It depends on the average distances amogns the points</param>
        /// <param name="nlag">Number of Lags. Typically is 12.</param>
        /// <param name="gdist">Get distaces in an array</param>
        /// <param name="gsemi">Get semivariances in an array </param>
        public void CalculateLags(
            double lag, int nlag, out ArrayList gdist, out ArrayList gsemi)
        {
            ////0.min 1.Max 2. Dist sum, 3. Semi sum, 4 number of points
            this.mlags = new double[nlag, 5];
            gdist = new ArrayList();
            gsemi = new ArrayList();
            double auxR = 0;
            for (int i = 0; i < nlag; i++)
            {
                this.mlags[i, 0] = auxR;
                this.mlags[i, 1] = auxR + lag;
                this.mlags[i, 2] = 0;
                this.mlags[i, 3] = 0;
                this.mlags[i, 4] = 0;
                auxR += lag;
            }

            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = i; j < this.npoints; j++)
                {
                    if (i != j)
                    {
                        double tdist = 0;
                        double tsemi = 0;
                        tdist = this.Mdistances(i, j);
                        tsemi = this.MsemiVarEmp(i, j);
                        for (int k = 0; k < nlag; k++)
                        {
                            if ((tdist >= this.mlags[k, 0]) && (tdist < this.mlags[k, 1]))
                            {
                                this.mlags[k, 2] += tdist;
                                this.mlags[k, 3] += tsemi;
                                this.mlags[k, 4] += 1;
                            }
                        }
                    }
                }
            }

            for (int i = 0; i < nlag; i++)
            {
                if (this.mlags[i, 4] > 0)
                {
                    this.mlags[i, 2] = this.mlags[i, 2] / this.mlags[i, 4];
                    this.mlags[i, 3] = this.mlags[i, 3] / this.mlags[i, 4];
                    gdist.Add(this.mlags[i, 2]);
                    gsemi.Add(this.mlags[i, 3]);
                }
                else
                {
                    this.mlags[i, 2] = this.mlags[i, 0] + (0.5 * lag);
                }
            }
        }

        /// <summary>
        /// Generate the arrays of distances and semivariances for drawing the theoretical semivariance
        /// </summary>
        /// <param name="gdistF">Array distance</param>
        /// <param name="gsemiT">Array semivariances</param>
        public void CalculateModelGraphic(
            out ArrayList gdistF, out ArrayList gsemiT)
        {
            double div = this.maux / 30;
            gdistF = new ArrayList();
            gsemiT = new ArrayList();
            for (int i = 0; i < 30; i++)
            {
                gdistF.Add(div * i);
                gsemiT.Add(this.mmodel.GetValue(div * i));
            }
        }



        /// <summary>
        /// This method generates an ArraList with distances for drawing the semivariogram
        /// </summary>
        /// <returns> An array with the distances</returns>
        public ArrayList Distances()
        {
            ArrayList mdist = new ArrayList();
            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = i; j < this.npoints; j++)
                {
                    if (j != i)
                    {
                        mdist.Add(this.Mdistances(i, j));
                    }
                }
            }

            return mdist;
        }

        /// <summary>
        /// This method generates an ArraList with semivariances for drawing the semivariogram
        /// </summary>
        /// <returns>An array with all semivariances</returns>
        public ArrayList SemiVariances()
        {
            ArrayList msemi = new ArrayList();
            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = i; j < this.npoints; j++)
                {
                    if (j != i)
                    {
                        msemi.Add(this.MsemiVarEmp(i, j));
                    }
                }
            }

            return msemi;
        }

        /// <summary>
        /// This method gets the value depending on the number of points to search nPn, and the position of the point
        /// </summary>
        /// <param name="nPn">Number of points to execute the interpolation</param>
        /// <param name="x">Coordinate X of the point to calculate</param>
        /// <param name="y">Coordinate Y of the point to calculate</param>
        /// <returns>Value found after interpolation</returns>
        public override double[] Interpolate(
            double x, double y, double z, bool include)
        {
            List<Kpoint> cpointsF;

            this.Search.ListPointAvailables(out cpointsF, x, y, z, include);

            if (cpointsF.Count == 1)
            {
                return new double[2]{cpointsF[0].W,double.NaN};
            }

            if (cpointsF.Count == 0)
            {
                return new double[2] { double.NaN, double.NaN };
            }


            int nPn = cpointsF.Count;
            double value = 0;
            vectorG = new double[nPn + 1];
            weights = new double[nPn + 1];
            vectorG[nPn] = 1.0;


            double[,] mdistF = new double[nPn + 1, nPn + 1];
            double[,] msemEF = new double[nPn + 1, nPn + 1];
            double[,] msemTF = new double[nPn + 1, nPn + 1];
            double[,] inverMatrixF = new double[nPn + 1, nPn + 1];


            this.CalculateDistancesAround(ref cpointsF, ref mdistF, ref msemEF);

            try
            {
                this.GenerateMatrixProximity(nPn, ref msemTF, ref inverMatrixF, ref mdistF);
            }
            catch
            {
                return new double[2] { double.NaN, double.NaN };
            }


            for (int i = 0; i < nPn; i++)
                vectorG[i] = this.mmodel.GetValue(Utils.GetDistance(x, y, z,cpointsF[i]));


            for (int i = 0; i < nPn + 1; i++)
            {
                weights[i] = 0;
                for (int j = 0; j < nPn + 1; j++)
                {
                    weights[i] += inverMatrixF[i, j] * vectorG[j];
                }
            }

            for (int i = 0; i < nPn; i++)
            {
                value += weights[i] * cpointsF[i].W;
            }


            double error = 0;

            for (int i = 0; i <= nPn; i++)
            {
                error += weights[i] * vectorG[i];
            }

          //  if (error > 0)
                error = Math.Sqrt(error);


            return new double[2]{value,error};
        }




        /// <summary>
        /// This method calculates the distance between points that are close to a point
        /// </summary>
        /// <param name="npoints">Number of points around</param>
        /// <param name="cpoints">Matrix with coordinates</param>
        /// <param name="mdistances">Distances between points</param>
        /// <param name="msemiVarEmp">Semivariance of point selected</param>
        private void CalculateDistancesAround(
             ref List<Kpoint> cpoints, ref double[,] mdistances, ref double[,] msemiVarEmp)
        {
            for (int i = 0; i < cpoints.Count; i++)
            {
                for (int j = 0; j < cpoints.Count; j++)
                {
                    mdistances[i, j] = Utils.GetDistance(cpoints[i], cpoints[j]);

                    msemiVarEmp[i, j] = Utils.GetSemivariance(cpoints[i], cpoints[j]);

                    if (mdistances[i, j] > this.maux)
                    {
                        this.maux = mdistances[i, j];
                    }

                    if (msemiVarEmp[i, j] > this.mauxV)
                    {
                        this.mauxV = msemiVarEmp[i, j];
                    }
                }
            }
        }




        /// <summary>
        /// This method calculate the distances between all points
        /// </summary>
        /// 
        private void CalculateDistances()
        {
            double mdis = 0d;
            double msem = 0d;
            Parallel.For(0, this.npoints, i =>
             {
                 Parallel.For(0, this.npoints, j =>
                 {
                     mdis = this.Mdistances(i, j);
                     msem = this.MsemiVarEmp(i, j);

                     if (mdis > this.maux)
                     {
                         this.maux = mdis;
                     }

                     if (msem > this.mauxV)
                     {
                         this.mauxV = msem;
                     }
                 });
             });
        }

        /// <summary>
        /// Generate initial matrix for interpolation using the point around of unknown value
        /// </summary>
        /// <param name="nPn">Number of point used</param>
        /// <param name="msemiVarTeo">Matrix semivariances</param>
        /// <param name="inva">Matrix inverse</param>
        /// <param name="mdistances">Matrix of distaces</param>
        private void GenerateMatrixProximity(
            int nPn, ref double[,] msemiVarTeo, ref double[,] inva, ref double[,] mdistances)
        {
            for (int i = 0; i <= nPn; i++)
            {
                msemiVarTeo[i, nPn] = 1.0;
                inva[i, nPn] = 1.0;
            }

            for (int i = 0; i <= nPn; i++)
            {
                msemiVarTeo[nPn, i] = 1.0;
                inva[nPn, i] = 1.0;
            }

            msemiVarTeo[nPn, nPn] = 0;
            inva[nPn, nPn] = 0;
            for (int i = 0; i < nPn; i++)
            {
                for (int j = 0; j < nPn; j++)
                {
                    msemiVarTeo[i, j] = this.mmodel.GetValue(mdistances[i, j]);
                    inva[i, j] = this.mmodel.GetValue(mdistances[i, j]);// msemiVarTeo[i, j];
                }
            }

            //Matriz matriz = new Matriz(inva);
            Statistics.Matrix matrix = Statistics.Matrix.FromDoubleArray(inva).Invert();

            inva = Statistics.Matrix.ToDoubleArray(matrix);


          
           
            ////int info = 0;
            //// matinv.matinvreport rep = new matinv.matinvreport();
            //// matinv.rmatrixinverse(ref inva, nPn + 1, ref info, ref rep);
        }






        public override double[] CrossValidation()
        {
            double[] output = new double[5];
            List<Kpoint> cpointsCross = new List<Kpoint>();
            double SSerror = 0;
            double sumAllValues = 0;
            double quad = 0; ;
            double sumError = 0;
            int numberUsed = 0;
            double averageError=0;
            double meanStdPredError=0;
            double quadMeanStdPredError = 0;

            for (int drop = 0; drop < this.NumPoints; drop++)
            {
                cpointsCross = new List<Kpoint>();
                for (int i = 0; i < this.NumPoints; i++)
                    if (i != drop)
                        cpointsCross.Add(this.Kpoints[i]);

             
                Kpoint value = this.Kpoints[drop];
                double[] inter = this.Interpolate(value.X, value.Y,value.Z, false);
                double error = inter[1];
                if (
                    !inter[0].Equals(double.NaN) && !inter[1].Equals(double.NaN))
                {
                    value.Zinte = inter[0];
                    value.Error = inter[1];

                    if (value.Error != 0)
                    {

                    value.SdtError = value.Zinte - value.W;

                    averageError += value.Error;

                    sumError += value.SdtError;
                    meanStdPredError += value.SdtError / value.Error;
                    quadMeanStdPredError += ((value.SdtError / value.Error) * (value.SdtError / value.Error));
                
                    sumAllValues += value.Zinte;
                    quad += (value.Zinte * value.Zinte);
                    SSerror += (value.SdtError * value.SdtError);
                    numberUsed++;
                    }
                }
            }

            if (numberUsed>0)
            {
            output[0] = Math.Sqrt(SSerror / numberUsed); //RMS
            output[1] = sumError / numberUsed; //Mean Prediction errors
            output[2] = Math.Sqrt(averageError / numberUsed); //average krigin std error
            output[3] = meanStdPredError / numberUsed; // Meand sdt prediction error
            output[4] = Math.Sqrt(quadMeanStdPredError / numberUsed); //RMS sdt prediction error
}
            else
            {
              return null;
            }
            return output;
        }


        public IRaster SurfaceError(DotSpatial.Data.Extent inExtents, double cellSize, DotSpatial.Projections.ProjectionInfo projection, double noDataValue, string outGridPath)
        {
            IRaster raster;
            // this.np=this.Search.MaxNumPoints;
            CreateGridFromExtents(inExtents, cellSize, projection, noDataValue, outGridPath, out raster);

            int nr = raster.NumRows;
            int nc = raster.NumColumns;

            for (int row = 0; row < nr; row++)
            {
                for (int col = 0; col < nc; col++)
                {
                    Coordinate cor = this.GetXYCell(raster, col, row);
                    double val = this.Interpolate(cor.X, cor.Y,cor.Z, true)[1];
                    raster.Value[row, col] = val;
                }
            }
            raster.Save();
            raster.Close();
            return raster;
        }

        public IRaster SurfaceError(DotSpatial.Data.IRaster mask, DotSpatial.Data.Extent inExtents, double cellSize, DotSpatial.Projections.ProjectionInfo projection, double noDataValue, string outGridPath)
        {
            IRaster raster;
            // this.np=this.Search.MaxNumPoints;
            CreateGridFromExtents(inExtents, cellSize, projection, noDataValue, outGridPath, out raster);

            int nr = raster.NumRows;
            int nc = raster.NumColumns;
            
            for (int row = 0; row < nr; row++)
            {
                for (int col = 0; col < nc; col++)
                {
                    if (mask.Value[row, col] == 1)
                    {
                        Coordinate cor = this.GetXYCell(raster, col, row);
                        double val = this.Interpolate(cor.X, cor.Y,cor.Z, true)[1];
                        raster.Value[row, col] = val;
                    }
                    else
                        raster.Value[row, col] = noDataValue;
                }
            }
            raster.Save();
            raster.Close();
            return raster;
        }

        public DataSourceInterpolators GetDataSource()
        {
            return (DataSourceInterpolators)this;
        }


        public void SetKrigingModel(IKrigingModel model)
        {
            this.KrigingModel = model;
        }
    }
}
